MK440: Project 2 - Kiwi Bubbles

# Load data and packages
library("dummies")

## dummies-1.5.6 provided by Decision Patterns

library("AER")

## Loading required package: car

## Loading required package: carData

## Loading required package: lmtest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: sandwich

## Loading required package: survival

library("plotly")

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library('RColorBrewer')

library("data.table")
library("mlogit")

## Loading required package: Formula

library("gmnl")

## Loading required package: maxLik

## Loading required package: miscTools

## 
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
## 
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/

set.seed(0)
setwd("~/desktop/Pricing/Project2")
data=fread("kiwi_bubbles_P2.csv")
data=data[!(data$price.KB==99),]
data=data[!(data$price.KR==99),]
data=data[!(data$price.MB==99),]

Load packages and data. Do some data cleaning by removing data that contains prices at 99 (pre-processing).

# Logit model with regression(No segmentation)
 ## Get average price
avg_priceKB=mean(data$price.KB)
avg_priceKR=mean(data$price.KR)
avg_priceMB=mean(data$price.MB)
table(avg_priceKB,avg_priceKR,avg_priceMB)

## , , avg_priceMB = 1.34558500323206
## 
##                  avg_priceKR
## avg_priceKB       1.37808661926309
##   1.3813316095669                1

 ## Define demand function
demand=function(priceKB,priceKR,priceMB,coef){
  probKB=exp(coef[1]+coef[4]*priceKB)/(1+exp(coef[1]+coef[4]*priceKB)+exp(coef[2]+coef[4]*priceKR)+exp(coef[3]+coef[4]*priceMB))
  probKR=exp(coef[2]+coef[4]*priceKR)/(1+exp(coef[1]+coef[4]*priceKB)+exp(coef[2]+coef[4]*priceKR)+exp(coef[3]+coef[4]*priceMB))
  probMB=exp(coef[3]+coef[4]*priceMB)/(1+exp(coef[1]+coef[4]*priceKB)+exp(coef[2]+coef[4]*priceKR)+exp(coef[3]+coef[4]*priceMB))
  return(cbind(probKB,probKR,probMB))
}
  ## Define profit function
uc=0.5 #unit cost
profit=function(priceKB,priceKR,priceMB,coef){
  profitKB=demand(priceKB, priceKR,priceMB, coef)[,1]*(priceKB-0.5)*1000
  profitKR=demand(priceKB, priceKR,priceMB, coef)[,2]*(priceKR-0.5)*1000
  return(cbind(profitKB,profitKR))
}

 ## Run regression
mlogitdata=mlogit.data(data,id="id",varying=4:7,choice="choice",shape="wide")
mle= gmnl(choice ~  price, data = mlogitdata)
summary(mle)

## 
## Model estimated on: Thu Feb 13 02:53:03 2020 
## 
## Call:
## gmnl(formula = choice ~ price, data = mlogitdata, method = "nr")
## 
## Frequencies of categories:
## 
##       0      KB      KR      MB 
## 0.41564 0.18035 0.20039 0.20362 
## 
## The estimation took: 0h:0m:0s 
## 
## Coefficients:
##                Estimate Std. Error z-value  Pr(>|z|)    
## KB:(intercept)  4.25316    0.32821  12.959 < 2.2e-16 ***
## KR:(intercept)  4.36240    0.32945  13.241 < 2.2e-16 ***
## MB:(intercept)  4.20440    0.31331  13.419 < 2.2e-16 ***
## price          -3.73793    0.23671 -15.791 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Optimization of log-likelihood by Newton-Raphson maximisation
## Log Likelihood: -1909
## Number of observations: 1547
## Number of iterations: 4
## Exit of MLE: gradient close to zero

coef=mle$coefficients

 ## Probability
prob=demand(avg_priceKB,avg_priceKR,avg_priceMB,coef)
prob

##                   probKB    probKR    probMB
## KB:(intercept) 0.1753666 0.1979974 0.1908972

 ## Elasticity
prob=demand(avg_priceKB,avg_priceKR,avg_priceMB,coef)
probKB=prob[1]
probKR=prob[2]
probMB=prob[3]

KBelas=-coef[4]*avg_priceKB*(1-probKB)
KRelas=-coef[4]*avg_priceKR*(1-probKR)
MBelas=-coef[4]*avg_priceMB*(1-probMB)

KBKRcrosselas=-coef[4]*avg_priceKR*probKR
KBMBcrosselas=-coef[4]*avg_priceMB*probMB
KRKBcrosselas=-coef[4]*avg_priceKB*probKB
KRMBcrosselas=-coef[4]*avg_priceMB*probMB
MBKBcrosselas=-coef[4]*avg_priceKB*probKB
MBKRcrosselas=-coef[4]*avg_priceKR*probKR

elasticity=data.frame('KB'=c(0,0,0),'KR'=c(0,0,0),'MB'=c(0,0,0))
rownames(elasticity)=c('KB','KR','MB')
elasticity[1,1]=KBelas
elasticity[1,2]=KBKRcrosselas
elasticity[1,3]=KBMBcrosselas
elasticity[2,1]=KRKBcrosselas
elasticity[2,2]=KRelas
elasticity[2,3]=KRMBcrosselas
elasticity[3,1]=MBKBcrosselas
elasticity[3,2]=MBKRcrosselas
elasticity[3,3]=MBelas
elasticity

##           KB       KR        MB
## KB 4.2578474 1.019923 0.9601564
## KR 0.9054743 4.131270 0.9601564
## MB 0.9054743 1.019923 4.0695469

 ## Calculate optimal price
uc=0.5#Unit cost
priceMB=1.43
aux=seq(0,3,0.01)
pricespace=expand.grid(aux,aux)
profitmat=matrix(0L,nrow(pricespace),1)

for (i in 1:nrow(pricespace)){
  profitmat[i]=sum(profit(pricespace[i,1],pricespace[i,2],priceMB,coef))  
}

xaxis=list(title="P^{KB}")
yaxis=list(autorange = "reversed",title="P^{KR}")
zaxis=list(title="Profit")
p=plot_ly(x=pricespace[,1],y=pricespace[,2],z=as.numeric(profitmat),
          type="scatter3d",mode="markers",
          marker = list(color = as.numeric(profitmat), colorscale = c('#FFE1A1', '#683531'), showscale = TRUE))%>%
  layout(scene=list(xaxis=xaxis,yaxis=yaxis,zaxis=zaxis))%>%
  config(mathjax = 'cdn')
p

max(profitmat)

## [1] 393.4082

 ## Optimal price for KB and KR
OptimalPrice=pricespace[which(profitmat==max(profitmat)),]
OptimalPrice#KB=1.16 KR=1.16

##       Var1 Var2
## 35033 1.16 1.16

After calculating the elasticity and cross elatsticity, we found that a product’s own elasticity is the same as its cross-elasticity with other products.These results make sense. We can find out why in its definition. A product j’s own elasticity equals to −β1Pj(1 − Pr(y = j | P)) while a product j’s cross elasticity is defined as −β1Pj’Pr(y = j’| P) where j’ is j’s rival product. Under this scenario, it means that, for instance, KR’s cross elasticity with KB is −β1PKB(1 − Pr(y = KB | P)) because KB can be treated as its rival, which is the same as KB’s own elasticity,−β1PKB(1 − Pr(y = KB | P)). It is easy to understand, that with rival product j’ pricing up, product j would sell more.

# Q4 Logit model with segmentation
 ## Load packages and data
data=fread("kiwi_bubbles_P2.csv")
data=data[!(data$price.KB==99),]
data=data[!(data$price.KR==99),]
data=data[!(data$price.MB==99),]
demo = fread("demo_P2.csv",stringsAsFactors = F)
 ## Use Kmeans to segment
N = 329
demo_cluster = kmeans(x=demo[, 2:18], centers = 7, nstart = 1000)

 ##Combine cluster identity into the raw data
cluster_id = data.frame(id = demo$id)
cluster_id$cluster = demo_cluster$cluster
data = merge(data, cluster_id, by = "id", all.x = T)
data$cluster[is.na(data$cluster)] = 8
seg.share = c( table(demo_cluster$cluster),N - sum(table(demo_cluster$cluster))) / N 
seg.share

##          1          2          3          4          5          6          7 
## 0.09422492 0.14285714 0.12462006 0.16109422 0.11550152 0.11854103 0.10334347 
##            
## 0.13981763

coef.est = data.frame(segment = 1:8, intercept.KB = NA, intercept.KR = NA, 
                      intercept.MB = NA, price.coef = NA) 
#Write a for-loop. 
for (seg in 1:8) {
  # During each loop, pick subset of data of consumers from each segment.
  data.sub = subset(data, cluster == seg)
  
  #Using that data, the rest remains the same.
  mlogitdata=mlogit.data(data.sub,id="id",varying=4:7,choice="choice",shape="wide")
  
  #Run MLE.
  mle= gmnl(choice ~  price, data = mlogitdata)
  mle
  #Store the outcome in the coef.est matrix.
  coef.est[seg, 2:5] = mle$coefficients
}

 ##Define demand function
demandKB=function(priceKB,priceKR,priceMB,para){
  prob=exp(para[1]+para[4]*priceKB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
  return(prob)
}
demandKR=function(priceKB,priceKR,priceMB,para){
  prob=exp(para[2]+para[4]*priceKR)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
  return(prob)
}
demandMB=function(priceKB,priceKR,priceMB,para){
  prob=exp(para[3]+para[4]*priceMB)/(1+exp(para[1]+para[4]*priceKB)+exp(para[2]+para[4]*priceKR)+exp(para[3]+para[4]*priceMB))
  return(prob)
}
demandKR_noKB=function(priceKR,priceMB,para){
  prob=exp(para[1]+para[3]*priceKR)/(1+exp(para[1]+para[3]*priceKR)+exp(para[2]+para[3]*priceMB))
  return(prob)
}
 ##Define aggregate demand and profit function
agg_choice=function(priceKB,priceKR,priceMB) {
  
  agg_choice_KB = seg.share[1]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[1,2:5]))+
    seg.share[2]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[2,2:5]))+
    seg.share[3]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[3,2:5]))+
    seg.share[4]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[4,2:5]))+
    seg.share[5]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[5,2:5]))+
    seg.share[6]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[6,2:5]))+
    seg.share[7]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[7,2:5]))+
    seg.share[8]*demandKB(priceKB,priceKR,priceMB,as.numeric(coef.est[8,2:5]))
    
  agg_choice_KR = seg.share[1]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[1,2:5]))+
    seg.share[2]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[2,2:5]))+
    seg.share[3]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[3,2:5]))+
    seg.share[4]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[4,2:5]))+
    seg.share[5]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[5,2:5]))+
    seg.share[6]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[6,2:5]))+
    seg.share[7]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[7,2:5]))+
    seg.share[8]*demandKR(priceKB,priceKR,priceMB,as.numeric(coef.est[8,2:5]))
    
  return(cbind(agg_choice_KB,agg_choice_KR))
}

agg_choice1=function(priceKR,priceMB) {
  
  agg_choice_KR = seg.share[1]*demandKR_noKB(priceKR,priceMB,as.numeric(coef.est[1,3:5]))+
    seg.share[2]*demandKR_noKB(priceKR,priceMB,as.numeric(coef.est[2,3:5]))+
    seg.share[3]*demandKR_noKB(priceKR,priceMB,as.numeric(coef.est[3,3:5]))+
    seg.share[4]*demandKR_noKB(priceKR,priceMB,as.numeric(coef.est[4,3:5]))+
    seg.share[5]*demandKR_noKB(priceKR,priceMB,as.numeric(coef.est[5,3:5]))+
    seg.share[6]*demandKR_noKB(priceKR,priceMB,as.numeric(coef.est[6,3:5]))+
    seg.share[7]*demandKR_noKB(priceKR,priceMB,as.numeric(coef.est[7,3:5]))+
    seg.share[8]*demandKR_noKB(priceKR,priceMB,as.numeric(coef.est[8,3:5]))
  return(agg_choice_KR)
}

profit_Seg_noKB_F=function(priceKR,priceMB){
  profitKR=agg_choice1(priceKR,priceMB)*(priceKR-0.5)*1000
  return(profitKR)
}

pricespace=seq(0,3,0.01)
ProfitSeg_noKB = profit_Seg_noKB_F(pricespace,1.43)
KRMaxPnoKBSeg = pricespace[ProfitSeg_noKB==max(ProfitSeg_noKB)]
profit_Seg_noKB_F(1.06,1.43)

##        1 
## 294.2091

profit_Seg1_wKB=function(priceKB,priceKR,priceMB){
  profitKB=agg_choice(priceKB,priceKR,priceMB)[,1]*(priceKB-0.5)*1000
  profitKR=agg_choice(priceKB,priceKR,priceMB)[,2]*(priceKR-0.5)*1000
  return(cbind(profitKB,profitKR))
}
 ##Find optimal price
aux=seq(1,3,0.01)
pricespace=expand.grid(aux,aux)
profitmatSeg=matrix(0L,nrow(pricespace),1)

for (i in 1:nrow(pricespace)){
  profitmatSeg[i]=sum(profit_Seg1_wKB(pricespace[i,1],pricespace[i,2],1.43))  
}

pricespace[profitmatSeg==max(profitmatSeg)]

## [1] 1.15 1.18

sum(profit_Seg1_wKB(1.15,1.18,1.43))

## [1] 394.0433

 ##Elasticity
Nominator_KB=0
Nominator_KR=0
Nominator_MB=0
for (i in 1:8) {
  KB=seg.share[i]*demandKB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*coef.est[i,5]*(1-demandKB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5])))
  Nominator_KB=KB+Nominator_KB
}
for (i in 1:8) {
  KR = seg.share[i]*demandKR(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*coef.est[i,5]*(1-demandKR(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5])))
  Nominator_KR=KR+Nominator_KR
}
for (i in 1:8) {
  MB = seg.share[i]*demandMB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*coef.est[i,5]*(1-demandMB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5])))
  Nominator_MB=MB+Nominator_MB
}

Denominator_KB=0
Denominator_KR=0
Denominator_MB=0
for (i in 1:8) {
  KB=demandKB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*seg.share[i]
  Denominator_KB=KB+Denominator_KB
}
for (i in 1:8) {
  KR=demandKR(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*seg.share[i]
  Denominator_KR=KR+Denominator_KR
}
for (i in 1:8) {
  MB=demandMB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*seg.share[i]
  Denominator_MB=MB+Denominator_MB
}

KB_Seg_Elas=-avg_priceKB*Nominator_KB/Denominator_KB
KR_Seg_Elas=-avg_priceKR*Nominator_KR/Denominator_KR
MB_Seg_Elas=-avg_priceMB*Nominator_MB/Denominator_MB

 ##Cross-elasticity
Nominator_KBKR=0
Nominator_KBMB=0
Nominator_KRKB=0
Nominator_KRMB=0
Nominator_MBKB=0
Nominator_MBKR=0

for (i in 1:8) {
  KBKR=seg.share[i]*demandKB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*coef.est[i,5]*demandKR(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))
  Nominator_KBKR=KBKR+Nominator_KBKR
}
for (i in 1:8) {
  KBMB=seg.share[i]*demandKB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*coef.est[i,5]*demandMB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))
  Nominator_KBMB=KBMB+Nominator_KBMB
}
for (i in 1:8) {
  KRKB=seg.share[i]*demandKR(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*coef.est[i,5]*demandKB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))
  Nominator_KRKB=KRKB+Nominator_KRKB
}
for (i in 1:8) {
  KRMB=seg.share[i]*demandKR(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*coef.est[i,5]*demandMB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))
  Nominator_KRMB=KRMB+Nominator_KRMB
}
for (i in 1:8) {
  MBKB=seg.share[i]*demandMB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*coef.est[i,5]*demandKB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))
  Nominator_MBKB=MBKB+Nominator_MBKB
}
for (i in 1:8) {
  MBKR=seg.share[i]*demandMB(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))*coef.est[i,5]*demandKR(avg_priceKB,avg_priceKR,avg_priceMB,as.numeric(coef.est[i,2:5]))
  Nominator_MBKR=MBKR+Nominator_MBKR
}

KBKR_Seg_Cross_Elas=-avg_priceKR*Nominator_KBKR/Denominator_KB
KBMB_Seg_Cross_Elas=-avg_priceMB*Nominator_KBMB/Denominator_KB
KRKB_Seg_Cross_Elas=-avg_priceKB*Nominator_KRKB/Denominator_KR
KRMB_Seg_Cross_Elas=-avg_priceMB*Nominator_KRMB/Denominator_KR
MBKB_Seg_Cross_Elas=-avg_priceKB*Nominator_MBKB/Denominator_MB
MBKR_Seg_Cross_Elas=-avg_priceKR*Nominator_MBKR/Denominator_MB

Seg_elasticity=data.frame('KB'=c(0,0,0),'KR'=c(0,0,0),'MB'=c(0,0,0))
rownames(Seg_elasticity)=c('KB','KR','MB')
Seg_elasticity[1,1]=KB_Seg_Elas
Seg_elasticity[1,2]=KBKR_Seg_Cross_Elas
Seg_elasticity[1,3]=KBMB_Seg_Cross_Elas
Seg_elasticity[2,1]=KRKB_Seg_Cross_Elas
Seg_elasticity[2,2]=KR_Seg_Elas
Seg_elasticity[2,3]=KRMB_Seg_Cross_Elas
Seg_elasticity[3,1]=MBKB_Seg_Cross_Elas
Seg_elasticity[3,2]=MBKR_Seg_Cross_Elas
Seg_elasticity[3,3]=MB_Seg_Elas
Seg_elasticity

##           KB        KR        MB
## KB 4.4745579 0.9859642 1.0915897
## KR 0.8324154 3.8223969 0.8835077
## MB 0.9853480 0.9446296 4.2396141

We have decided to put the dataset into 8 segments before the number of each segment become too small. After we calculated the own elasticities and cross elasticities for the segmentation case, we have observed lower own-elasticities for KB(3.028), KR(2.96) and a higher own-elasticities for MB(4.95) than the non-segmentation case. We also found the cross elasticity between KB with MB, KR with MB much lower. We found the cross elasticity between KB and KR(1.24) are higher than KB and MB(0.50) meaning KB and KR are close substitute and KB and MB are not close substitute. Looking at the beta zero of each segments. We have found segment 2,3,8 has higher preference on KB and for the case of Kiwi Regular without Kiwi Bubbles and Mango Bubbles priced at 1.43, we will set optimal price for Kiwi Regular at 1.06 and the max profit is 294.20.

When we launch Kiwi Bubbles, the optimal price is 1.15 for Kiwi Regular and 1.18 for Kiwi Bubbles. The profit of Kiwi has increased from 294.2091 to 394.0433 after we launch Kiwi Bubbles. By launching KB, we are able to charge a higher price for KR and we are taking market share from MB so that we are earning a higher profit and lowering our competitors profit.

Moreover, the own elasticity of KB and KR when there is no segmentation are 4.25 and 4.13, after we do segmentation, they became 3.02 and 2.96. It means that customers are less sensitive to our price change, or in other words, our customers loyalty has been increased. On The contrary, competitor’s own elasticity has increased from 4.07 to 4.95, which means they are losing customers’ brand loyalty.

Thus, the model justified the launching of KB is correct.

# Q5 Understanding strategic responses
 ## Re-define demand and profit function
demand=function(priceKB,priceKR,priceMB,coef){
  probKB=exp(coef[1]+coef[4]*priceKB)/(1+exp(coef[1]+coef[4]*priceKB)+exp(coef[2]+coef[4]*priceKR)+exp(coef[3]+coef[4]*priceMB))
  probKR=exp(coef[2]+coef[4]*priceKR)/(1+exp(coef[1]+coef[4]*priceKB)+exp(coef[2]+coef[4]*priceKR)+exp(coef[3]+coef[4]*priceMB))
  probMB=exp(coef[3]+coef[4]*priceMB)/(1+exp(coef[1]+coef[4]*priceKB)+exp(coef[2]+coef[4]*priceKR)+exp(coef[3]+coef[4]*priceMB))
  return(cbind(probKB,probKR,probMB))
}

profit=function(priceKB,priceKR,priceMB,coef){
  profitKB=demand(priceKB, priceKR,priceMB, coef)[,1]*(priceKB-0.5)*1000
  profitKR=demand(priceKB, priceKR,priceMB, coef)[,2]*(priceKR-0.5)*1000
  profitMB=demand(priceKB, priceKR,priceMB, coef)[,3]*(priceMB-0.5)*1000
  return(cbind(profitKB,profitKR,profitMB))
}
 ## First round
aux=seq(0,3,0.01)
profitmat=c()
for (i in 1:length(aux)){
  profitmat[i]=profit(1.16,1.16,aux[i],coef)[,3]
}
max(profitmat)#177.6243

## [1] 177.6243

MB1=aux[which(profitmat==max(profitmat))]
MB1#=0.95

## [1] 0.95

aux=seq(0,3,0.01)
pricespace=expand.grid(aux,aux)
profitmat=matrix(0L,nrow(pricespace),1)

for (i in 1:nrow(pricespace)){
  profitmat[i]=sum(profit(pricespace[i,1],pricespace[i,2],0.95,coef)[,1:2])  
}
max(profitmat)#275.7644

## [1] 275.7644

KBKR1=pricespace[which(profitmat==max(profitmat)),]
KBKR1#1.04 1.04

##       Var1 Var2
## 31409 1.04 1.04

 ## Second Round
aux=seq(0,3,0.01)
profitmat=c()
for (i in 1:length(aux)){
  profitmat[i]=profit(1.04,1.04,aux[i],coef)[,3]
}
max(profitmat)#145.6836

## [1] 145.6836

MB2=aux[which(profitmat==max(profitmat))]
MB2#=0.91

## [1] 0.91

aux=seq(0,3,0.01)
pricespace=expand.grid(aux,aux)
profitmat=matrix(0L,nrow(pricespace),1)
for (i in 1:nrow(pricespace)){
  profitmat[i]=sum(profit(pricespace[i,1],pricespace[i,2],0.91,coef)[,1:2])  
}
max(profitmat)#262.2523

## [1] 262.2523

KBKR2=pricespace[which(profitmat==max(profitmat)),]
KBKR2#KB=1.03 KR=1.03

##       Var1 Var2
## 31107 1.03 1.03

 ## Third Round
aux=seq(0,3,0.001)
profitmat=c()
for (i in 1:length(aux)){
  profitmat[i]=profit(1.03,1.03,aux[i],coef)[,3]
}
max(profitmat)#143.0388

## [1] 143.0389

MB3=aux[which(profitmat==max(profitmat))]
MB3#=0.91

## [1] 0.911

aux=seq(0,3,0.01)
pricespace=expand.grid(aux,aux)
profitmat=matrix(0L,nrow(pricespace),1)
for (i in 1:nrow(pricespace)){
  profitmat[i]=sum(profit(pricespace[i,1],pricespace[i,2],0.91,coef)[,1:2])  
}
max(profitmat)#262.2523

## [1] 262.2523

KBKR3=pricespace[which(profitmat==max(profitmat)),]
KBKR3#KB=1.03 KR=1.03

##       Var1 Var2
## 31107 1.03 1.03

After two rounds’ price war, KB KR and MB converge. KB and KR’s prices would be 1.03 and MB’s price would be 0.91. This kind of actions are more likely to happen in real world because a rival would know when I price down to gain more profits and of course, rival would react to kind of situation. In the previous scenario, we assume that rival would not react but in Q5, the process is dynamic until KB KR and MB converge at a new equalibrium point. Clearly, this scenario conforms more to the economic principles.